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Abstract 

We present a novel Discontinuous Galerkin Finite Element Method for wave propagation problems. The method 
employs space-time Trefftz-type basis functions that satisfy the underlying partial differential equations and the re- 
spective interface boundary conditions exactly in an element-wise fashion. The basis functions can be of arbitrary high 
order, and we demonstrate spectral convergence in the L2-norm. In this context, spectral convergence is obtained with 
respect to the approximation error in the entire space-time domain of interest, i.e. in space and time simultaneously. 
Formulating the approximation in terms of a space-time Trefftz basis makes high order time integration an inherent 
property of the method and clearly sets it apart from methods, that employ a high order approximation in space only. 

Keywords: Discontinuous Galerkin Method, Finite Element Method, Trefftz Method, Higher order time integration, 
Electrodynamics, Wave propagation, 



1. Introduction 

In the context of Finite Element Methods (FEM) there exist two main ways of improving the numerical accu- 
racy. First, a computational mesh on which a problem is approximated can be refined (/z-refinement). Secondly, 
the order p of the approximating polynomials can be increased (p-refinement). Both strategies, their combinations 
(/jp-refinement) and their connection to the level of regularity of the solutions have been extensively studied (see 
e.g.EEl). 

Even higher accuracy can be obtained by applying problem-specific high order basis functions. We pursue this 
approach in the framework of the Discontinuous Galerkin Finite Element Methods (DG-FEM) ||3][4l[5][6), as it allows 
an almost unrestricted choice of bases. To date, there exist only few works on DG and DG-type methods using a 
problem specific basis in the frequency-domain, in particular the Ultra Weak Variational Formulation (UWVF) ||71E[9] 
and the Plane Wave Discontinuous Galerkin Method (PWDG) flOl fTTl which has been shown to be a generalization 
of the former one. In the time-domain, we are aware of lfl2l only. However, this approach differs significantly from 
our work as will be detailed below. Most commonly however, generic polynomials are employed, which does not lead 
to the optimal approximation accuracy for a given polynomial order and mesh size. 

In this article we present a highly accurate type of DG-FEM for time-domain applications. A distinguishing new 
attribute of the method is the use of Trefftz-type basis functions lfT3l[T4l[T5l . which, by definition, exactly satisfy the 
underlying partial differential equations and the respective interface boundary conditions in an element-wise fashion. 
The method is, hence, a Discontinuous Galerkin Trefftz Finite Element Method (DGT-FEM). 

We consider space- and time-dependent equations. As Trefftz functions are required to solve the equations exactly, 
at least in a local sense, they necessarily have to depend on time as well. Consequently, the DGT-FEM has to be 
formulated as a space-time FEM I l6l[T7ll . In a space-time method there is no formal distinction between the spatial 
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dimensions and the temporal one. Finite element spaces are given as (d + l)-dimensional functional spaces, where d 
denotes the number of spatial dimensions. 

Working in a space-time setting has two immediate implications: first, high order time integration is obtained 
in a straightforward and consistent manner by increasing the order of the approximating polynomials, and second, 
convergence rates describe the approximation error not only in space but in space-time. It is shown below that 
spectral convergence of the global approximation error in the space-time domain of interest is obtained under p- 
enrichment for sufficiently smooth data. This is in contrast to DG methods based on a discretization of space, where 
spectral convergence is attained with respect to the spatial discretization error only. We also show that the presence 
of inhomogeneous elements (i.e elements filled with different materials) does not affect the convergence rate if an 
appropriate local Trefftz basis is used. 

The remainder of this article is organized as follows. In Section|2]we derive the weak space-time DG formulation 
of time-dependent electromagnetic waves. The derivation is restricted to the (1 + l)-dimensional case. A generalization 
to higher spatial dimensions is outlined below and will be implemented in the near future. Section [3] introduces the 
space-time Trefftz basis functions that by necessity have to be vector-valued, with the electric and magnetic field 
as components. Trefftz bases are derived for homogeneously filled grid cells as well as for partially filled cells 
containing two different materials. Results obtained from numerical experiments are presented in Section|4] First, we 
demonstrate the physical correctness of the results. Then we investigate the convergence rates of the method both in 
space and time. Finally, we show that the method outperforms well-established methods such as the Finite Difference 
Time Domain Method (FDTD) and the DG method with leapfrog time stepping (DGL) in terms of efficiency. Section|5] 
concludes the paper and provides an outlook to the ongoing work of extending the method to higher dimensions in 
space. 



2. Time-Dependent One-Dimensional Waves 

This section contains a derivation of a weak space-time DG formulation of Maxwell's equations using vectorial 
basis functions. The resulting equations describe the propagation of electromagnetic waves in one spatial dimension. 
The procedure developed for electromagnetic waves in this section is applicable to waves of any nature. However, we 
shall focus on electromagnetic analysis for definiteness. 

2.1. Differential Formulation of the Problem in Continuum 

For a wave propagating in a given direction x with one-component fields E = E y and H = H z , the system of 
Maxwell's equations simplifies to 

d x E + d t {jiH) = and d x H + 8,{eE) = J, (1) 

where only the time-dependent equations are considered. Here, e and p. are given material parameters that can in 
general be functions of both position and time (in electrodynamics, these are the dielectric permittivity and the mag- 
netic permeability, respectively). J is a given source that can also depend on space x and time t. The electric field 
E e H l (Q) and the magnetic field H e H l (Q) are defined in the Sobolev space 

H\Q.) = {ue L 2 (Q) : Vm e L 2 (0)}, 

with L2(Q) being the space of square integrable functions over the space-time domain of interest 

CI = I, X I v c M 2 with I f = [f min ; f max ] c R and I v = [^ min ; * max ] c E. 

For highlighting the generality of the following derivation, we rewrite equations ([1} in a coordinate-independent 
divergence form 



d,Y I fi \ I E 
d r ' 1 I \ H 



o. 

(2) 



2 



By defining the material and differential operators 

*=( £ ?)>^(? o) andv =(i )' 

we can further contract equations (|2]) into the more compact form 

^■^■¥ = and V T ■ 77, • F = J. (3) 

Here F = (E, H) T e H 1 (£2) x H 1 (£2) is the field vector defined in the vectorial Sobolev space. Equations ((3} are subject 
to initial as well as boundary conditions. 

2.2. The Weak Form 

In order to obtain a weak formulation of (|3), we first cast it into the more abstract form 

£F = J, (4) 

where we used the source vector J = (0, J) T e LiiQ) x Lqiil) and the differential Maxwell operator _£ : //'(£!) x 
H l (Q.) — > L2(Cl) x L2(Q) defined as £, = (V T ?/ e , V 1 ?^). We then produce a weak formulation, by composing the inner 
product of both sides of (Jij) with a test function v = (v £ , v H ) T e //'(Q) x H l (Q.) as 

[XF,T] n = [J,v] n , (5) 

where [■, ]n is the inner product in LaiCl) x L2(^) defined as 

[F,y] £1 = (E,v E ) n + (H,v H ) n , 

and (•, )q is the standard L2 inner product in £1 Casting the weak formulation |5]) back into a more explicit form 
yields 

(V^F, y £ ) n + (V T ^F, v w ) n = (7, y H ) n . (6) 

2.3. The DG Formulation 

To obtain a DG formulation of equation |6]l the expression is transformed using integration by parts. This yields 

f v £ (7, e F)-ndO+ f y // (77 /J F)-ndQ-(V T -v £ ,7 7e F) n -(V T -v H ,/7 /J F) n = (y,y H ) n , (7) 
Jan Jon 

where n = (n,, n V ) T is the unit outward normal on the space-time domain boundary dQ.. More specifically, we consider 
a partition Q/, of the space-time domain of interest O on a Cartesian grid (not necessarily uniform) into K — N ■ M 
cells. Here denotes the total number of cells in the temporal direction whereas M denotes the total number of cells 
in the spatial direction. Each cell with index k = (n,m) is denoted by Q.^ where n is the temporal index and m the 
spatial index. The components of the test functions v E and v H are defined cell wise in their discretized Sobolev spaces 

H\Q. h ) = W e L 2 (Q) : VO t e Q. h , v\ Qt e ff 1 ^*)] with i = E,H. 

The vectorial test function is therefore defined in v € //'(O/,) x //'(Q/,). 

To generate a time-stepping scheme on £2/,, we will consider only two spatial layers of grid cells, defined at two 
consecutive time slabs labeled n and n + 1 (see. Fig. [TJ. It is inherent to our method that space and time inside one 
time slab are treated on an equal footing in the discretization procedure. We will apply (j7]) to a single grid cell £2 m ,„+i 
at time-level n + 1 rather than n, as our focus is on implicit schemes. The field vector F is separately approximated 
within each cell using a linear combination of basis functions u, viz. 

F„, m (*, t) = V f,(„X hm (x, t), (8) 
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Figure 1: Space-time cells for two subsequent time layers. 

where each basis function is vectorial, with the two components u£„, = (Un'm,Un'm) T 6 H l (Q., hm ) x //'(□„„,). Note 
that the approximation is of order p in space and time while applying the same number of Degrees of Freedom (DoF) 
as a corresponding approximation of space only. 

In contrast to most approaches, we are interested in employing Trefftz bases rather than the standard piecewise- 
polynomial approximations. However, in this part of the derivation we keep a general form of basis functions and 
implement the Trefftz basis in Section [3] 

Once a set of basis functions over the cell is chosen, the integral over C2„+i,m directly leads to the entries of the 
stiffness matrices 

S%£T )= f {U E ' p {d x u^) + u H ^{d x u^)\ dxdt, 

J J V 'n+l,m 

S™iT )= f [(eu^'id^+^idru"-")) dxdt. 

J J V 'n+l,m 

On the element level, the treatment of the integrals over the element boundary F„ + i m is more nuanced than in the 
global problem (|7b. These integrals represent surface fluxes that in the continuous problem must be the same for any 
two adjacent grid cells. In DG methods basis and test functions are defined with a strictly element-wise compact 
support. Therefore, the approximation is continuous within each element but discontinuous across any cell interfaces. 
At the interfaces the global approximation is double-valued at every point and interface fluxes from adjacent grid cells 
can differ. Flux continuity can be numerically imposed via a penalty term or, as is commonly done in DG methods, 
by combining weighted numerical approximations of each flux in the two adjacent cells. More specifically, we are 
using flux expressions that are centered in space and upwind in time. Let us write out all the fluxes in terms of basis 
and test functions of maximum order P m and Q m , respectively, in a given space-time cell Q n>m . We denote evaluations 
of the basis functions on the left and right hand cell boundary using left and right arrows (*—, — >), respectively. In 
temporal direction we indicate evaluations of the basis functions at the upper (i.e. towards the next time level) and 
lower boundary using up and down pointing arrows (T,D- The element boundary F„ + i m is dissected into the four 
element edges. The space flux in cell Q.„+\, m is composed from three terms. Cell £2„ +1 m itself contributes 

t q,P (space) _ £ f ( E,p,->H,q,^ H,p,-.£^,-t\ , _ 1 f ( ,, E 'P^ ,, H '1'^ + „ H ' P '^ u E ' q ^\ dt 
^n+\,m,m ~ 2 I \ m m m m / J 2 I \ / i 

In addition, the neighboring cell on the right Q„ + i „ 1+ i contributes with 

L? P ( S paoe) = 1 [ / E,p,<- H,q,-> H,pj- E q ,-,\ rf 
n+l,m,— > 2 J \ / +i 
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and to the neighboring cell on the left Q n+ i m _i with 

L ^( 8 pace) = _ 1 ft U E,P^ U H,^ H,p,—> E,q,*—\ rf 
«+l,m,«— 2 J \ — / +1 

The time flux of cell Q„ + i, m has only two contributions from cells Q n +i sW and Q nm 

L , P (space) = r l E ^ E,^ Hp* H qt \ 

T ?.P(time) _ _ f (E.P.I E,q,i H,p,\ H,q\\ , 

L «+l,m,| - J \ €U " U n+l + f XU » U n+1 " X " 

We can now rewrite Maxwell's equations in a fully discretized matrix form, which reads 

yq,p (time) (space) _ g?,p(time) _ g^,/? (space) \ w> (space) rp yq,p (space) „p _ j^,p(time) 

We further combine the stiffness and flux matrix expressions on the left and right hand sides into two global matrices 
A and B and formally obtain the implicit time update equation 

A/„ +1 = B/„ (9) 

for progressing from time layer n to a subsequent time layer n + l. 



3. Approximation with Trefftz Functions 



Trefftz basis functions are a set of problem-specific functions, which satisfy the underlying partial differential 
equations exactly within a grid cell. As they include features of the exact solution, Trefftz methods often yield 
very accurate numerical solutions. In this section we construct a Trefftz basis for electromagnetic wave propagation 
problems as stated in Section [2] For some scenarios, such as the propagation of a plane wave in a vacuum, Trefftz 
basis functions are relatively easy to construct |[T8l[T9l . However, we also address the case of partially filled elements, 
where the derivation of locally exact functions is more involved. 

3.1. Trefftz Functions for (1+1 )-Dimensional Maxwell's Equations 

In the following we apply Maxwell's equations ([TJ in natural units, with the permittivity and permeability of 
vacuum set to one. As a consequence, the speed of light in vacuum, v , is one as well. We omit cell indices (i.e. k, n 
and m) for the sake of simplicity from now on. 

Traveling waves of the form ip(x, t) = f(x + vf), where / is an arbitrary sufficiently smooth function and v the 
local speed of light, are well known to be solutions of the wave equation. Motivated by this, it is reasonable to include 
transport polynomials of the general form 



,H,p,± 



' ±{x + V t) 



(10) 



in the basis. The first and second component, u E and u H , of any basis function in ( 10 1 are representing the electric 
and magnetic field, respectively. The material parameters Z 



jfi/e and v = 1/ yfjli (i.e. intrinsic impedance and 
local speed of light respectively) enter the basis functions directly. For each order p, there are two waves in the basis: 
one moving leftward (with opposite signs of E and H) and the other one rightward (equal sign of E and H). We note 
that using separate scalar bases for approximating the electric and magnetic field individually does not yield a Trefftz 
method. Any individual basis function or a combination of them has to solve the governing equations. In the two field 
formulation that we apply, this is not guaranteed with separate scalar basis functions as can easily be seen by inserting 
transport polynomials of order r and s with r + s into ([TJ. The approximation of the field vector F is obtained as 



-Z Z 



7?,sign ..p,sign 



sign=+,— 



)■ 



(ii) 
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Figure 2: Transport polynomials of orders /) = 0, 1, 2, 3 in a vacuum (i.e. e = 1 and fi = 1) to describe waves traveling rightwards. 



In Figs. [2] and [3] we plot the basis functions for two different media in the master element. In both figures the cell 
size in spatial and temporal direction is two, however, while e = fi = 1 representing a vacuum is used in Fig. [2] Fig. [3] 
represents the situation for a dielectric with e = 4 and fi — I. In the first case, the direction of transport is along the 
cell diagonal, whereas in the second case, the direction of transport is rotated away from the diagonal to an angle of 
22.5 degrees. 





Figure 3: Transport polynomials of orders p = 0, 1, 2, 3 in a medium with 6 = 4 and fi = 1; corresponding to waves traveling rightwards. 

As mentioned above, Peterson et al. fT2 1 deal with a DG time-domain method based on transport polynomials as 
well. However, they discretize the second order wave equation, resulting in a weak formulation where both trial and 
test functions occur in the form of derivatives only. Therefore, all constant terms vanish rendering the stiffness matrix 
singular and requiring the application of regularization techniques. A reconstruction procedure that ensures point wise 
inter element continuity is used to define the element average, which is unnatural in the context of DG methods and 
makes the implementation of material interfaces much harder, if possible at all in higher dimensions, inter-element 
continuity is weakly imposed by means of Lagrange multipliers instead of numerical fluxes. Apart from being an 
unusual approach within a DG method, the multipliers are additional degrees of freedom, which have to be solved for. 
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3.2. Treffiz Basis Functions for Partially Filled Cells 

In the following, a Trefftz basis is derived for partially filled cells. To this end, we consider a space-time cell 
[-A x /2, Ajt/2] x [-A,/2, A,/2] with the origin placed at the center. Let there be a material interface at some point 
x = xo inside this cell; the medium velocities and intrinsic impedances on the left and the right side of the interface 
being v\ t2 and Zi j2 , respectively. We look for a Trefftz function corresponding, physically, to a traveling incident wave 
('inc') that is coming from the left and being reflected ('ref') / transmitted ('trs') at the interface. The fields in the 
considered cell are 

£inc = <Pinc(x-Vit); H inc = ^-tfi inc (x - Vit) , 

£ref = 0>ref(* + Vif); #ref = — <Pref (X + V t f) , 

Zl 

Eas = <Ptxs(x-V 2 t); //trs = 7T<Ptrs(x-V 2 t) . 

At the interface, these fields must be continuous 

Vine 0(1 - ViO - (firtf (XQ + ViO = (p as (x - V 2 t) , 



TT (fine (XO ~ V\t) + lfi ref (X + V\t)) = ^-(p m (*0 - V 2 t) ■ 

Zi z 2 
Solving for the reflected and transmitted waveforms, we obtain 

Zi -z 2 



^ref (XQ + ViO = Rfinc (xq ~ V\ t) with R = 

ftis (xo - v 2 f) = Tip inc (xq - v\f) with T = 



z, +z 2 ' 

Z 2 
Zi + z 2 ' 



The reflection and transmission coefficients R and T are thus expressed in the same way as in the frequency domain, 
which should be expected for frequency-independent material parameters. What remains is to find these waveforms 
at an arbitrary position x rather than just at xq. We therefore express the outgoing waveforms in terms of the incoming 
waveform which yields 

<p Ttf (X + ViO = </>ref (XQ + Vit + (X- X )) = <fi Ttf (x + v/) 

= R(p inc (x - v^') = Rtp inc (2x -x-v\t), 



for the reflected waves and 



<Ars ( x ~ v 2t) = (firs (X(l - V 2 t + (x - X )) = ^trs (x - V 2 t") 

= Tip inc (x - v^") = T(p inc {jl - ^jx + ^x- vifj, 



for the transmitted waves. Here the time variables of the reflected and transmitted waveforms contain an advance or 
delay respectively 

f , = f+ (x-xo) ^ f „ = f _(x-xo)_ 
v\ v 2 

A set of Trefftz functions can now be obtained for any order p by setting Ei nc (ip) = ip p and H mc (<p) = tp p /Z. An 
additional completely similar set can be generated for incident waves impinging on the interface boundary from the 
right. 
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4. Results 



In this section we show numerical results obtained for different scenarios with the DGT-FEM as introduced above. 
These include simulations of wave propagation in a vacuum, and two scenarios with a material interface. In the first 
case the interface is located in between two elements and in the second one inside an element. The accuracy of the 
method is subsequently investigated, and we compare it to well-established methods. 

4.1. Propagation of a Gaussian Wave in Vacuum 

Gaussian wave propagation provides a simple test scenario for DGT-FEM. As the waveform has high regularity, 
strictly increasing accuracy is expected for high approximation orders. In Fig. [4] we show the simulation of the 
propagation of a Gaussian shaped wave (initially heading leftwards). In this test case, we assume that the wave travels 
in a vacuum (i.e. e = fi = 1 and v = 1). The cell size is normalized to unity (i.e. A v = A, = 1) for convenience. At the 
spatial boundaries of the global space-time domain, perfect electric conductor (PEC) boundary conditions are applied. 
Fig. [4] shows the electric field of the Gaussian wave heading leftwards at an angle of 45 degrees in the space-time 
plane corresponding to a velocity v = 1. At the boundary the wave is reflected under a sign change of the electric field 
while maintaining its space-time angle of 45 degrees, thus recovering the expected physical behavior. 




Figure 4: Electric field of a one-dimensional Gaussian wave in a vacuum, simulated with DGT-FEM. The solution in the whole space-time domain 
of interest (x, t) e [-20, 20] X [0, 60] is displayed. 



4.2. Propagation of a Gaussian Wave in Cell-Wise Homogeneous Media 

We now modify the setup of the preceding section by filling the left part of the domain with a dielectric material 
with e — 4. The medium interface is placed between two cells at xq = -10. Fig. [5] illustrates the result in the 
space-time domain of interest. The wave is initialized in a vacuum and heads leftwards with a space-time angle of 45 
degrees towards the dielectric material. At time t = 20 the wave reaches the interface and is split into a reflected and 
a transmitted wave. The reflected wave heads back rightwards preserving its 45 degrees angle in order to be reflected 
at the right domain boundary. The numerical value of the amplitude of this wave is very close to the analytical value 
R = 1/3. For the transmitted wave the space-time angle changes to 22.5 degrees, which corresponds to a velocity of 
Vmed = 1/2; the amplitude of this wave is very close to the analytical value T = 2/3 again. Therefore, the expected 
reflection and transmission coefficients as well as wave velocities are recovered. For a wave impinging from medium 
to a vacuum we obtain R' = -1/3 and T = 4/3, which agrees with the expectations as well. 

4.3. Simulation of Inhomogeneously Filled Cells 

In Section [372] Trefftz basis functions for the case of material interfaces placed within a cell were derived. Fig. [6] 
shows the simulation of the propagation of a wave similar to the preceding example. While the media have exactly 
the same material properties as before, the medium interface is now placed at xq = -0.25; that is inside a cell. In this 
case, we also recover correct physical properties (i.e. reflection and transmission coefficients as well as local speed of 
light) as discussed before. 
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Figure 5: Electric field of a one-dimensional Gaussian wave with a medium interface at xq = -10, simulated with the DGT-FEM. The solution in 
the whole space-time domain of interest (x, t) e [-20, 20] X [0, 60] is displayed. 



We note that while the preceding examples can be handled in a very similar manner with non-Trefftz methods, 
this is not the case here. Typically material interfaces inside elements are either completely avoided, or the material 
properties are averaged in the respective element. More complex methods such as immersed boundaries or cut cells 
are in principle capable of handling this situation but they usually come at high numerical costs and exhibit sub- 
optimal performance. Using a Trefftz approach with the above derived basis, however, the wave behavior at the 
interior interface is modeled accurately. 




Figure 6: Electric field of a one-dimensional Gaussian wave with a medium interface at xq = -0.25, simulated with the DGT-FEM. The material 
interface crosses an element. The solution in the whole space-time domain of interest (x, t) e [-5, 5] X [0, 20] is displayed. 



4.4. Approximation Error and Efficiency 

In Section[T]we mentioned the advantage of Trefftz methods regarding accuracy compared to those using generic 
polynomials. While the results presented so far are of a more qualitative nature, we perform a quantitative error 
analysis in this section. To this end we first investigate the numerical error occurring during the propagation of a 
Gaussian wave over a large distance in a vacuum. 

In Fig. [7] the relative L2-error obtained with the DGT method (green markers) is compared to the error obtained 
using a centered DG discretization of space combined with a leapfrog time stepping scheme (DGL) (6] |20| (blue 
markers). As can be seen from the graph, the accuracy of the DGL scheme is limited by the second order time 
integrator. However, the new DGT-FEM achieves spectral convergence as is inferred from the straight line in the 
semi-logarithmic plot. Note that the errors indicated in the plot include spatial as well as temporal errors. This 
behavior is a remarkable result of DGT-FEM. For most methods spectral convergence is only achieved in the space 
domain, which imposes a strong limitation on the total error. To yield a better convergence of the global error higher 
order time stepping schemes (e.g. high order Runge-Kutta schemes) are usually employed, which increase the overall 
computational costs. In DGT-FEM this is achieved at no additional computational costs. In Fig. [8] we plot the relative 
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Figure 7: Lo norm of the relative error for the propagation of a Gaussian wave with DGT and DGL for different polynomial orders. 



^2 -projection-error for one cell in DGT rilled with two different materials. Spectral convergence is also achieved in 
this case as can be inferred from the straight line in the semi-logarithmic plot again. Therefore, spectral convergence 
even holds for partially filled cells in DGT. 
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Figure 8: Li norm of the relative projection error for basis functions of one DGT cell filled with two different materials as a function of the 
polynomial order. 

Next, we will separately verify the spatial and temporal convergence orders for a fixed approximation order p. 
To this end, we perform series of simulations where one of the grid parameters, Ax or Af is successively divided 
in halves while the other respective parameter is fixed. This way the convergence rates in space and time can be 
evaluated individually. The results are depicted in Fig. [9] and Fig.[lO] The plots show that the expected convergence 
rates (i.e. p + 1 for basis functions of order p) are obtained for cell size variation of both x and t. 



In Fig. 1 1 we plot the relative L2 error for three methods including DGT (green curves) as a function of computing 
time. The other two methods are the DGL method mentioned above (blue curves) and a Finite Difference time domain 
(FDTD) method (red curve) respectively. The setup is identical for all methods, and we only show the computing time 
consumed by the update procedure. The error behavior of spatial DG methods (such as the DGL method used here) 
has been discussed in great detail in ET1I221I . One directly observes that DGT becomes the most efficient method for 
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Figure 9: L2 norm of the relative error as obtained with DGT for different polynomial orders as a function of the spatial cell size Ax. 
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Figure 10: L2 norm of the relative error as obtained with DGT for different polynomial orders as a function of the temporal cell size At. 



orders p > 1. 



5. Summary and Outlook 

We have developed a theoretical framework for the description of (electromagnetic) wave propagation problems in 
the time-domain in terms of a space-time DG method. In this framework we have subsequently implemented problem- 
specific Trefftz-type basis functions for cells with homogeneous as well piecewise homogeneous media. High-order 
time integration is an inherent property of the method, and we demonstrated spectral convergence in space-time for 
wave propagation problems. As the Trefftz functions for inhomogeneous media reproduce the exact physical behavior 
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Figure 1 1 : Z/> norm of the relative error for the propagation of a Gaussian wave with DGT and DGL with different polynomial orders as well as 
FDTD as function of CPU time. 



at media interfaces, spectral convergence is obtained even in the presence of partially filled elements. A performance 
investigation showed that the method outperforms established time-domain methods such as FDTD and DGL. 

The method is applicable to higher dimensional systems (i.e. (2+1) and (3+1)), and we are currently in the process 
of extending our implementation to problems governed by Maxwell's equations 



V x E + d,(jM) = and V x H - d t (eE) = J. 



(12) 



The fundamental difference to the (l+l)-dimensional case is that the directions of propagation are not known a 
priori (i.e. there are infinitely many), as has been the case in the (l+l)-dimensional case (i.e. waves were only able to 
propagate leftwards, and rightwards). Instead, we approximate the fields locally, within each homogeneous element, 
by plane waves of the form 



( n(k r - k \t) F 



kxfl (k-r - k-\t) p ) 



v z 



(13) 



where k is the direction of transport unit vector and n the unit vector characterizing the polarization of the respective 
basis function. Approximation of fields by plane waves is a well known and common technique. Its accuracy in the 
frequency domain has been studied rigorously by several research groups @[TTJ[23][24] but the mathematical analysis 
is very intricate and complex. We are not aware of similar mathematical results in the time domain and will therefore 
apply ( |T~3] > heuristically but with a high level of confidence in its validity, supported by the numerical calculations 
below. It can be expected that the accuracy of plane-wave approximation for higher dimensional systems will depend 
on the order of the polynomial basis, on the number of directions of propagation employed and, to a lesser extent, on 
the choice of these directions (see 1241 1. More generally, approximation properties of Trefftz bases ("T-sets") have 
been studied by Herrera and others (e.g. [25, 26 1). 

To get some quantitative estimates, we approximated a wave of the form ^(x, y, t) = ((jc— 1) cos(a') + y sin(o') - 1) 1 
with a = n/5 using (\3\ with increasing order and an increasing number of directions. The result shown in Fig. [T2] 
verifies that the approximation error depends on both parameters. In the case of five directions and order seven, the 
wave direction and order is matched, and the error approaches zero. Here, we benefit from the works on UWVF and 
PWDG J8] |9] [H] [lOll regarding the choice of directions. 
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Figure 12: L2 norm of the relative error of a (2+l)-dimensional seventh order initialization as function of the number of directions. 
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